﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;


namespace bstest
{
    class cal
    {
		//计算主程序
        public static void calmain(int calMode, dataIO.ImpellerData[] data,int number, double length, double elasticity, double density)
        {
            double[] shape = new double[51];
            double[] ans = new double[100];
            double frequency=0;
            int n = 50,count=0;
            elasticity = elasticity * Math.Pow(10, 6);
            density = density / Math.Pow(10, 6);
            switch (calMode)//根据计算类型，分配不同的计算类
            {
                case 1: calModelIterative(data, number, length, elasticity, density, n,out shape,out frequency); break;
                case 2: calRayleigh(data, number, length, elasticity, density, n, out shape, out frequency); break;
                case 3: calTransferMatrix(data, number, length, elasticity, density,n,out count,out ans); break;
            }
            if (calMode < 3) dataIO.WriteAns(n, shape, frequency, length);
            else dataIO.WriteAnsMat(count,ans);
        }
        
		//积分法的拟合
        public static void fit(int n, dataIO.ImpellerData[] data,double length,out double[,] a,int mode)
        {
            int m = 0, i, j;
            a = new double[15, 51];
            for (i = 0; i <= n; i++) for (j = 1; j < 15; j++) a[j,i] = 0; 
            double nowlength = 0;
            double dx = length / n,xa,xj;
            for (i = 0; i <= n; i++)
            {
                if (nowlength > data[m].l * 1.001) m++;
                xa = data[m].a + (nowlength - data[m].l) / (data[m + 1].l - data[m].l) * (data[m + 1].a - data[m].a);
                xj = data[m].j + (nowlength - data[m].l) / (data[m + 1].l - data[m].l) * (data[m + 1].j - data[m].j);
                if (mode==1)
                {
                    a[1, i] = Math.Pow(nowlength / length, 2);
                    a[2, i] = xa;
                }
                else
                {
                    a[1, i] = xa * dx;
                    a[2, i] = dx;
                }
                a[8, i] = xj;
                nowlength += dx; 
            }
        }

		//牛顿-科特斯积分
        public static void NewtonCotes(ref double[,] a,int n, int x)
        {
            int i;
            for (i = 1; i <= n; i++) a[x - 1, i] = a[x - 2, i] + a[x - 2, i - 1];
            double temp = 0;
            if (x<8)
                for (i=n;i>=0;i--)
                {
                    a[x , i] = temp;
                    temp += a[x - 1, i];
                }
            else
                for (i=0;i<=n;i++)
                {
                    temp += a[x - 1, i];
                    a[x, i] = temp;
                }
        }

		//振型迭代法的计算
        public static void calModelIterative(dataIO.ImpellerData[] data, int number, double length, double elasticity, double density, int n,out double[] shape,out double frequency)
        {
            int i;
            double[,] a = new double[15, n+1];
            shape = new double[n+1];
            double amax,t,precision=0.001;
            bool f = true;
            fit(n,data,length,out a,1);
            while(f)
            {
                for (i = 0; i <= n; i++) a[3, i] = a[1, i] * a[2, i];
                NewtonCotes(ref a, n, 5);NewtonCotes(ref a, n, 7);
                for (i = 0; i <= n; i++) a[9, i] = a[7, i] / a[8, i];
                NewtonCotes(ref a, n, 11);NewtonCotes(ref a, n, 13);
                f = false;amax = 0;
                for (i = 0; i <= n; i++) if (a[13, i] > amax) amax = a[13, i];
                for (i=0;i<=n;i++)
                {
                    a[14, i] = a[13, i] / amax;
                    t = Math.Abs((a[14, i] - a[1, i]) / a[14, i]);
                    if (t > precision) f = true;
                    a[1, i] = a[14, i];
                }
            }
            for (i = 0; i <= n; i++) shape[i] = a[14,i];
            frequency = 4 / (Math.Pow(length / n, 2)) * Math.Sqrt(elasticity*1000/(density * a[13, n]))/(2*Math.PI);   
        }

		//雷利法的计算
        public static void calRayleigh(dataIO.ImpellerData[] data, int number, double length, double elasticity, double density, int n, out double[] shape, out double frequency)
        {
            int i;
            double[,] a = new double[15, n + 1];
            shape = new double[n + 1];
            double amax=0,sum1=0,sum2=0,g=9.8;
            fit(n, data, length, out a,2);
            for (i = 0; i <= n; i++) a[1, i] = a[1, i] * density;
            for (i = 0; i <= n; i++) a[3, i] = a[1, i] / a[2, i];
            NewtonCotes(ref a, n, 5); NewtonCotes(ref a, n, 7);
            for (i = 0; i <= n; i++) a[9, i] = a[7, i] / a[8, i];
            NewtonCotes(ref a, n, 11); NewtonCotes(ref a, n, 13);
            for (i = 0; i <= n; i++) if (a[13, i] > amax) amax = a[13, i];
            for (i = 0; i <= n; i++)
            {
                a[14, i] = g / elasticity * a[13, i] * Math.Pow(length / n, 4) / 16;
                shape[i] = a[13, i] / amax;
                sum1 += a[1,i] * a[14, i];
                sum2 += a[1, i] * Math.Pow(a[14,i],2);
            }
            frequency = Math.Sqrt(g * sum1 / sum2) / (2 * Math.PI) * Math.Pow(g,1.5);
        }

		//试凑法的矩阵拟合
        public static void fitforMat(int n, dataIO.ImpellerData[] data, double length, out dataIO.ImpellerData[] a)
        {
            int m = 0, i;
            double nowlength = 0;
            double dx = length / n;
            a = new dataIO.ImpellerData[n+1];
            for (i = 0; i <= n; i++)
            {
                a[i].l = dx;
                if (nowlength > data[m].l * 1.001) m++;
                a[i].a = data[m].a + (nowlength - data[m].l) / (data[m + 1].l - data[m].l) * (data[m + 1].a - data[m].a);
                a[i].j = data[m].j + (nowlength - data[m].l) / (data[m + 1].l - data[m].l) * (data[m + 1].j - data[m].j);
                nowlength += dx;
            }
        }

		//残矩的计算
        public static double calR(double omega,int n,double density, double elasticity,dataIO.ImpellerData[] a)
        {
            int i, j, k;
            double ej,x1,x2,R,m;
            double[,,] h = new double[5, 5, n + 1];
            double[] z1 = new double[5];double[] z2 = new double[5];
            for (i=0;i<=n;i++){
		        m=a[i].l* a[i].a* density;
                ej= elasticity * a[i].j;
		        h[1,1,i]=1;h[1,2,i]=a[i].l;h[1,3,i]=Math.Pow(a[i].l,2)/(2*ej);h[1,4,i]=h[1,3,i]*a[i].l/3;
		        h[2,1,i]=0;h[2,2,i]=1;h[2,3,i]=a[i].l/ej;h[2,4,i]=h[1,3,i];
		        h[3,1,i]=0;h[3,2,i]=0;h[3,3,i]=1;h[3,4,i]=a[i].l;
		        h[4,1,i]= Math.Pow(omega,2)*m;h[4,2,i]=h[4,1,i]*a[i].l;
		        h[4,3,i]=h[4,2,i]*a[i].l/(2*ej);h[4,4,i]=h[4,3,i]*a[i].l/3+1;
	        }
	        for (i=1;i<5;i++){z1[i]=0;z2[i]=0;}
            z1[3]=1;z2[4]=1;
	        for (i=1;i<=n;i++){
		        for (j=1;j<5;j++){
			        x1=0;x2=0;
			        for (k=1;k<5;k++){x1+=h[j,k,i]* z1[k];x2+=h[j,k,i]* z2[k];}
			        z1[j]=x1;z2[j]=x2;
		        }
	        }
	        R=z1[4]-z1[3]/z2[3]* z2[4];
	        return R;
        }

		//二分查找程序
        public static void find(double w1, double w2, double x1, double x2, ref double[] ans, ref int count, int n, double density, double elasticity, dataIO.ImpellerData[] a)
        {
            double mid = (w1 + w2) / 2,tht= 0.00001;
            double xmid = calR(mid, n, density, elasticity, a);
            if (Math.Abs(xmid) < tht) { ans[count] = mid; count++; }
            else
            {
                if (x1 * xmid < 0 && Math.Abs(xmid - x2) > tht) find(w1, mid, x1, xmid,ref ans,ref count, n, density, elasticity, a);
                if (x2 * xmid < 0 && Math.Abs(xmid - x1) > tht) find(mid, w2, xmid, x2, ref ans, ref count, n, density, elasticity, a);
            }
        }

		//传递矩阵法的计算
        public static void calTransferMatrix(dataIO.ImpellerData[] data, int number, double length, double elasticity, double density, int n,out int count,out double[] ans)
        {
            dataIO.ImpellerData[] a = new dataIO.ImpellerData[n + 1];
            fitforMat(n,data,length,out a);
            double step = 20,omega = step, x1, x2, tht = 0.00001;
            count = 0;
            ans = new double[100]; 
            while (omega < 10000)
            {
                x1 = calR(omega,n,density,elasticity,a); x2 = calR(omega + step,n,density,elasticity,a);
                if (Math.Abs(x1) < tht) { ans[count] = omega; count++; }
                else if (x1 * x2 < 0) find(omega, omega + step, x1, x2,ref ans,ref count, n, density, elasticity, a);
                omega += step;
            }
        }
    }
}
